個票データの集計処理例

関数定義ほか

単日データを処理するための関数を定義しておきます。

関数名	処理概要	必要なパッケージ
lagdiff(n)	単日データの前日差を求める	dplyr
ma7(n)	単日データの7日間（1週間）移動平均を求める	zoo
ma28(n)	単日データの28日間（4週間）移動平均を求める	zoo
ma(n, k)	単日データのk日間移動平均を求める	zoo

lagdiff <- function(n) {        # 前日差を求める関数
  n - dplyr::lag(n, default = 0L)
}

ma7 <- function(n) {            # 移動平均(7日)を求める関数
  zoo::rollmeanr(n, k = 7L, na.pad = TRUE)
}

ma28 <- function(n) {           # 移動平均(28日)を求める関数
  zoo::rollmeanr(n, k = 28L, na.pad = TRUE)
}

ma <- function(n, k = 7L) {     # 移動平均(k日)を求める関数
  zoo::rollmeanr(n, k = k, na.pad = TRUE)
}

累計は Base R の cumsum 関数を利用します。

　
その他、グラフで繰り返し使う文字列を変数として定義しておきます。

subtitle <- paste0("Generated @", lubridate::now())
caption <- "Data Source: covid19japan.com"

　

Import & Tidy

個票データの集計に限らず、データをインポートした際には目的に応じて各変量（変数）の変数型を変更します。特に水準ごとに層別処理を行いたい場合には因子型に変換しておくと便利です。また、結合したいデータと名称や体系を合わせておくこともポイントです。

　

都道府県データ

都道府県に関するデータは Gist で公開しているデータを用います。データ都道府県コード順に並んでいますので、この並びを利用して forcats パッケージを用いて因子化します。Base R の as.factor 関数を用いると水準がアルファニューメリック順に並べ変えられてしまう点に注意してください。

prefs <- "https://gist.githubusercontent.com/k-metrics/9f3fc18e042850ff24ad9676ac34764b/raw/f4ea87f429e1ca28627feff94b67c8b2432aee59/pref_utf8.csv" %>% 
  readr::read_csv() %>% 
  dplyr::mutate(
    # Googleの予測データと結合するためにコード体系を合わせる
    japan_prefecture_code = paste0("JP-", `コード`)
  ) %>% 
  dplyr::select(
    # Googleの予測データと結合するために名称を変更する
    japan_prefecture_code,　prefecture_name = pref,
    # 日本語の変数名は扱いにくいので英語名に変更する
    pref = `都道府県`, region = `八地方区分`, pops = `推計人口`
  ) %>% 
  dplyr::mutate(
    # 水準ごとに表示させるために因子化する（あらかじめデータをコード順に
    # 並べておくことが因子化の際のポイントのひとつ）
    japan_prefecture_code = forcats::fct_inorder(japan_prefecture_code),
    pref = forcats::fct_inorder(pref),
    region = forcats::fct_inorder(region),
    pops = as.integer(pops)
  )

prefs

　

Covid19Japan個票データ

個票データは Covid19Japan.com のデータリポジトリ で公開されているデータを用います。データは参照しているソースが更新されると自動的に更新されますので、取得タイミングによっては全てのデータが揃っていない可能性があります。
JSON形式データの各項目については データフォーマット を参照してください。

df <- "https://raw.githubusercontent.com/reustle/covid19japan-data/master/docs/patient_data/latest.json" %>% 
  # JSON形式をデータフレームとして読み込む
  jsonlite::fromJSON() %>% 
  dplyr::select(patientId, date = dateAnnounced, gender,
                pref = detectedPrefecture, patientStatus, knownCluster,
                confirmedPatient, ageBracket,
                deceasedDate, deceasedReportedDate) %>% 
  dplyr::filter(confirmedPatient == TRUE) %>% 
  dplyr::mutate(
    date = lubridate::as_date(date),
    gender = forcats::as_factor(gender),
    pref = stringr::str_to_lower(pref),
    patientStatus = forcats::as_factor(patientStatus),
    cluster = dplyr::if_else(!is.na(knownCluster), TRUE, FALSE),
    ageBracket = forcats::as_factor(ageBracket),
    deceasedDate = lubridate::as_date(deceasedDate),
    deceasedReportedDate = lubridate::as_date(deceasedReportedDate)
  ) %>% 
  # 都道府県データと結合
  dplyr::left_join(prefs, by = c("pref" = "prefecture_name")) %>% 
  dplyr::select(-pref) %>% 
  dplyr::rename(pref = pref.y) %>% 
  # 因子型の欠損値を水準化する
  dplyr::mutate(
    japan_prefecture_code = forcats::fct_explicit_na(japan_prefecture_code,
                                                     na_level = "JP-48"),
    pref = forcats::fct_explicit_na(pref, na_level = "空港検疫"),
    region = forcats::fct_explicit_na(region, na_level = "空港検疫"),
    gender = forcats::fct_explicit_na(gender, na_level = "非公表"),
    ageBracket = forcats::fct_explicit_na(ageBracket, na_level = "非公表"),
    patientStatus = forcats::fct_explicit_na(patientStatus,
                                             na_level = "Unknown")
  ) 

# %>% 
#   dplyr::filter(date < lubridate::today())

df

　
作成したデータフレームを確認します。意外と欠損値（n_missing項参照）が多いデータであることが分かります。

df %>% 
  skimr::skim()

Data summary

Name	Piped data
Number of rows	647000
Number of columns	14
_______________________
Column type frequency:
character	2
Date	3
factor	6
logical	2
numeric	1
________________________
Group variables	None

Variable type: character

skim_variable	n_missing	complete_rate	min	max	empty	n_unique	whitespace
patientId	0	1	1	16	0	647000	0
knownCluster	644495	0	3	88	0	233	0

Variable type: Date

skim_variable	n_missing	complete_rate	min	max	median	n_unique
date	0	1	2020-01-15	2021-05-10	2021-01-16	468
deceasedDate	646619	0	2020-02-13	2021-03-16	2020-05-08	152
deceasedReportedDate	646670	0	2020-02-13	2020-10-17	2020-05-16	131

Variable type: factor

skim_variable	n_missing	complete_rate	ordered	n_unique	top_counts
gender	0	1	FALSE	3	非公表: 537328, M: 61370, F: 48302
patientStatus	0	1	FALSE	9	Unk: 644465, Hos: 1261, Dec: 373, Hom: 315
ageBracket	0	1	FALSE	13	非公表: 537425, 20: 29433, 30: 19042, 40: 16089
japan_prefecture_code	0	1	FALSE	48	JP-: 147235, JP-: 90196, JP-: 56035, JP-: 40059
pref	0	1	FALSE	48	東京都: 147235, 大阪府: 90196, 神奈川: 56035, 埼玉県: 40059
region	0	1	FALSE	9	関東地: 298870, 近畿地: 157330, 中部地: 64081, 九州地: 57408

Variable type: logical

skim_variable	n_missing	complete_rate	mean	count
confirmedPatient	0	1	1	TRU: 647000
cluster	0	1	0	FAL: 644495, TRU: 2505

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
pops	2850	1	7601.01	4208.69	560	5107	7537	9177	13822	▇▅▆▇▇

　

Data Wrangling

　

単純集計

最初に単純集計を行ってみます。

　

全国集計

全国集計は日付データ（data）に基づく日次の単純集計です。厚生労働省のオープンデータにおける陽性者数データに相当します。単純集計は dplyr::group_by と dplyr::summrise を用いることで簡単に計数できます。ただし、データがない日は集計結果として現れませんので tidyr::complete で補完しておく必要があります。
その後、集計結果を元に前日差、累計、移動平均を計算します。

japan_daily <- df %>% 
  dplyr::group_by(date) %>% 
  dplyr::summarise(n = dplyr::n()) %>% 
  dplyr::ungroup() %>% 
  tidyr::complete(
    date = seq.Date(from = min(date), to = max(date), by = "day"),
    fill = list(n = 0L)
  ) %>% 
  dplyr::mutate(
    diff = lagdiff(n),   # 前日差
    cum = cumsum(n),     # 累計
    ma7 = ma7(n),        # 移動平均（7日）
    ma28 = ma28(n)       # 移動平均（28日）
  )

japan_daily

　

変数	変数の意味	変数型	備考
date	陽性判定者の報告日	Date	陽性判定日とは異なる
n	陽性判定者数	Integer
diff	同前日差	Integer	初日は陽性判定者数に同じ
cum	同累計	Integer
ma7	同移動平均（7日）	Double
ma28	同移動平均（28日）	Double	4週間（1ヶ月相当）

　

集計計算は dplyr::count を利用しても構いません。

df %>% 
  dplyr::count(date) # group_by(date) %>% summarize(n = n()) と等価

　

各変数別集計

dplyr::group_by で指定する変数を変えることで各変数の水準ごとの集計ができます。これは、前述の skimr::skim によるサマリーと同値になりますので下記以外の集計は省略します。

df %>% 
  dplyr::group_by(gender) %>% 
  dplyr::summarise(n = dplyr::n()) %>% 
  dplyr::ungroup()

　

クロス集計

クロス集計は日付と地方区分、日付と都道府県などのように複数の変数における水準ごとの集計です。tableなどを利用しても集計可能ですが、ggplot2 パッケージによる可視化を考慮し、ここではデータフレームを用いた集計を行います。
　
以下のような構造を取る個票データをクロス集計するとします。

df <- data.frame(
  date,   # Date型の日付データ
  key,    # グルーピングのためのデータ（例：都道府県、地方区分など）
  ...     # 順不同
)

基本的には単純集計と同じで dplyr::group_by と dplyr::summarize で集計します。補完も単純集計と同じく tidyr::complete で行うことができます。

x <- df %>% 
  dplyr::group_by(date, key) %>%      # クロス集計対象を指定する
  dplyr::summarise(n = dplyr::n()) %>%  # n()は個数をカウントする関数
  dplyr::ungroup() %>%                  # 最後にungroupするのがポイント
  tidyr::complete(                      # 暗黙の欠損を補完する
    date = seq.Date(from = min(date), to = max(date), by = "day"), key,
    fill = list(n = 0L)              # 個票がない=陽性者ゼロ
  )

　

クロス集計も dplyr::count で同様に処理できます。

df %>% 
  dplyr::count(date, gender) 

                     # group_by(date, gender) %>% summarize(n = n()) と等価

　

クロス集計後に前日差や累計などを算出するには tidyr::nest,tidyr::unnest と purrr パッケージを組み合わせた処理が確実かつ高速です。

x %>% 
  dplyr::group_by(key) %>%            # 上記の日次集計を元に各種計算を行う場合
  tidyr::nest() %>%                   # nest & unnest 処理を利用すると確実
  dplyr::mutate(                      # ネストしたデータはpurrrで処理する
    diff = purrr::map(data, ~ lagdiff(.$n)),   # 前日差
    cum = purrr::map(data, ~ cumsum(.$n)),     # 累計
    ma7 = purrr::map(data, ~ ma7(.$n)),        # 移動平均（7日）
    ma28 = purrr::map(data, ~ ma28(.$n))       # 移動平均（28日）
  ) %>% 
  tidyr::unnest(cols = c(data, diff, cum, ma7, ma28))

　
purrr パッケージの使い方がよく分からない場合は以下の方法でも可能です。

x %>% 
  dplyr::group_by(key) %>%            # ここがポイント
  dplyr::mutate(
    diff = lagdiff(n),
    cum = cumsum(n),
    ma7 = ma7(n),
    ma28 = ma28(n)
  ) %>% 
  dplyr::ungroup()

　

地方区分別集計

前述のクロス集計処理を daily_aggregate 関数として定義しているので、これを用いて集計を行っています。

region_daily <- df %>% 
  daily_aggregate(date, region)   # 前述の集計処理を関数化したもの

region_daily

　

都道府県別集計

地方区分集計と同様に集計します。

pref_daily <- df %>% 
  daily_aggregate(date, pref)

pref_daily

　

年代別集計

ageBracket_daily <- df %>% 
  daily_aggregate(date, ageBracket)

ageBracket_daily

　
年代別の地方区分別集計ならびに都道府県別集計は割愛します。理由はVisualize章のグラフを見て頂ければ分かると思います。

　

クラスター別集計

cluster_daily <- df %>% 
  daily_aggregate(date, cluster)

cluster_daily

　
クラスター別の地方区分別集計ならびに都道府県別集計は割愛します。理由はVisualize章のグラフを見て頂ければ分かると思います。

　

Visualize

作成した集計データを可視化します。 　

陽性判定者数

　

全国集計

japan_daily %>% 
  dplyr::filter(date == max(date))

単日と累計が二桁違うので、縦二軸のグラフを描きます。

title <- "【全国】陽性者数(単日)"
xlab <- ""
ylab <- ""
sec_scale <- 50   # 縦二軸のスケーリング

japan_daily %>% 
  ggplot2::ggplot(ggplot2::aes(x = date)) + 
    ggplot2::geom_bar(ggplot2::aes(y = n), stat = "identity", width = 1.0,
                      fill = "dark gray", alpha = 1.0) + 
    ggplot2::geom_line(ggplot2::aes(y = ma7), linetype = "solid",
                       colour = "gray10", size = 0.5) + 
    # 第二軸を利用するグラフを描画する際はスケーリング調整する必要あり
    ggplot2::geom_line(ggplot2::aes(y = cum / sec_scale),
                       colour = "dark green", size = 1.0) +
    # 横軸表示の指定
    ggplot2::scale_x_date(date_breaks = "1 month", date_labels = "%y/%m") + 
    # 二軸表示のための軸属性の指定
    ggplot2::scale_y_continuous(
      # 第一軸のラベル（スケールは自動調整）
      name = "陽性者数（灰）・移動平均（黒）",
      # 第二軸の指定（第一軸にあわせたスケーリング）
      sec.axis = ggplot2::sec_axis(~ . * sec_scale,
                                   name = "累積陽性者数（濃緑）")) + 
    # 縦軸の装飾
    ggplot2::theme(
      axis.text.y.left = ggplot2::element_text(colour = "gray10"),
      axis.line.y.left = ggplot2::element_line(colour = "gray10"),
      axis.text.y.right = ggplot2::element_text(colour = "dark green"),
      axis.line.y.right = ggplot2::element_line(colour = "dark green")
    ) +
    # グラフ全体のラベル指定
    ggplot2::labs(title = title, subtitle = subtitle, caption = caption,
                  x = xlab, y = ylab)

　
以降、描画コードは割愛しますが、このように個票を任意の変数の水準ごとに集計すると様々なデータの見かたができるようになります。

　

同前日差

　

同ヒートマップ

japan_daily %>% 
  dplyr::mutate(weekday = lubridate::wday(date, label = TRUE),
                weeks = lubridate::epiweek(date)) %>% 
  dplyr::mutate(weeks = dplyr::if_else(date > "2021-01-02",
                                       weeks + max(weeks), weeks)) %>% 
  # print() %>% 
  ggplot2::ggplot(ggplot2::aes(x = weekday, y = weeks, fill = n)) + 
    ggplot2::geom_tile() + 
    ggplot2::scale_fill_gradient(low = "white", high = "red") + 
    # ggplot2::coord_equal() + 
    # ggplot2::coord_flip() + 
    ggplot2::labs(title = "陽性者数ヒートマップ", x = "", y = "",
                  subtitle = subtitle, caption = caption)

　

同箱ひげ図

japan_daily %>% 
  dplyr::filter(n > 0) %>% 
  dplyr::mutate(weekday = lubridate::wday(date, label = TRUE)) %>% 
  ggplot2::ggplot(ggplot2::aes(x = weekday, y = n)) + 
    ggplot2::geom_violin(ggplot2::aes(fill = weekday), alpha = 0.25) +
    ggplot2::stat_boxplot(geom = "errorbar", width = 0.1) + 
    ggplot2::stat_summary(ggplot2::aes(group = 1, colour = weekday),
                          fun.y = median, geom = "line") + 
    ggplot2::geom_boxplot(width = 0.1) +
    ggplot2::labs(title = "曜日別陽性者数分布", x = "", y = "",
                  subtitle = subtitle, caption = caption)

　
　
参考までに箱ひげ図で示されている各分位数（陽性者ゼロの日を除き小数点以下は丸めてあります）は以下の通りです。

　

japan_daily %>% 
  dplyr::mutate(month = dplyr::if_else(lubridate::year(date) > 2020,
                                       lubridate::month(date) + 12,
                                       lubridate::month(date))) %>% 
  dplyr::filter(n > 0) %>% 
  dplyr::mutate(weekday = lubridate::wday(date, label = TRUE)) %>% 
  ggplot2::ggplot(ggplot2::aes(x = weekday, y = n)) + 
    ggplot2::geom_violin(ggplot2::aes(fill = weekday), alpha = 0.25) +
    ggplot2::stat_summary(ggplot2::aes(group = 1, colour = weekday),
                          fun.y = median, geom = "line") + 
    ggplot2::geom_boxplot(width = 0.1) +
    ggplot2::facet_wrap(~ month, nrow = 12, scales = "fixed", dir = "v") + 
    ggplot2::labs(title = "曜日別・月別陽性者数分布", x = "", y = "",
                  subtitle = subtitle, caption = caption)

　

地方区分別集計

　

　

同累計

　

同移動平均（7日）

　

同移動平均（28日）

　

同前日差

　

同ヒートマップ

　

同箱ひげ図（曜日別）

　

同箱ひげ図

subset <- region_daily %>% dplyr::mutate(key = region)
title <- "【地方別】陽性者数分布"
xlab <- ""
ylab <- ""
ncol <- 3

subset %>% 
  dplyr::filter(n > 0) %>% 
  ggplot2::ggplot(ggplot2::aes(x = reorder(key, n, median),
                               y = n, colour = key)) + 
    ggplot2::stat_boxplot(geom = "errorbar", width = 0.3) + 
    ggplot2::geom_boxplot() +
    ggplot2::theme(legend.position = 'none') + 
    ggplot2::coord_flip() + 
    ggplot2::labs(title = title, x = "", y = "",
                  subtitle = subtitle, caption = caption)

　

都道府県別

　

　

同前日差

　

同ヒートマップ

　

同箱ひげ図

　

同箱ひげ図

　

同箱ひげ図（直近3ヶ月）

　

年代別

　

クラスター別

　

対数化

地方区分（累計）

　

同移動平均（7日）

　

同移動平均（28日）

　

緊急事態宣言地域比較

週平均

ヒートマップ

　

---

CC 4.0 BY-NC-SA, Sampo Suzuki